/********************************************************************
Name: Figure 3: The Effect of Democratic House Member on
	Anti-Sanctuary Votes
Author: Dan Thompson
Date: August 2019
Notes:
********************************************************************/

* Set the working directory 
gl root = "$sanctuary_root"
gl root = "YOUR_WORKING_DIRECTORY"

* Bring in the house RD analysis data
use "$root/data/house immigration voting rd analysis data.dta", clear

* Set up the variables for storing the regression output, 
* predictions, and prediction error
gen yv = ""
gen b = .
gen se = .
gen n = .
gen y_hat = .
gen pred_error = .

* Loop over the roll call votes
levelsof year_vote, local(year_votes)
local i = 1
qui foreach yv in `year_votes' {
	
	* Store the vote number and year
	replace yv = "`yv'" if _n==`i'

	* Run the regression
	reg vote treat rv rv2 rv3 treat_rv treat_rv2 treat_rv3 ///
		if year_vote=="`yv'" & abs(rv)<25, robust

	* Store the key regression output
	replace b = _b[treat] if _n==`i'
	replace se = _se[treat] if _n==`i'
	replace n = e(N) if _n==`i'

	* Get the predicted values and error bars for the regression for plotting
	predict y_hat_`yv' if abs(rv)<25 & year_vote=="`yv'"
	predict pred_error_`yv' if abs(rv)<25 & year_vote=="`yv'", stdp
	replace y_hat = y_hat_`yv' if abs(rv)<25 & year_vote=="`yv'"
	replace pred_error = pred_error_`yv' if abs(rv)<25 & year_vote=="`yv'"
	drop y_hat_`yv' pred_error_`yv'
	
	* Increment the local indicating where I'm storing the regression output
	local ++i
}

* Construct the standard errors
gen pred_err_lower = y_hat - 1.96*pred_error
gen pred_err_upper = y_hat + 1.96*pred_error

* Make the bin means for the plot
gsort year_vote -treat rv
by year_vote: gen row = _n if treat==1
gsort year_vote treat -rv
by year_vote: replace row = _n if treat==0
gen bin = ceil(row/20)
egen mean_vote = mean(vote), by(year_vote bin treat)
egen mean_rv = mean(rv), by(year_vote bin treat)
gen mean_dem_vote_per = mean_rv+50

* Loop over the votes
levelsof year_vote, local(year_votes)
local i = 1
qui foreach yv in `year_votes' {
	
	* Limit the data to just one vote and within the bandwidth
	preserve
	keep if year_vote=="`yv'" & abs(rv)<25
	
	* Extract the year of the vote
	local year = substr("`yv'",1,4)
	
	* Run the RD plot
	twoway (scatter mean_vote mean_dem_vote_per, mc(gs7) msi(small)) ///
		(line pred_err_upper dem_vote_per if !treat, lc(gs10) lp(dash)) ///
		(line pred_err_upper dem_vote_per if treat, lc(gs10) lp(dash)) ///
		(line pred_err_lower dem_vote_per if !treat, lc(gs10) lp(dash)) ///
		(line pred_err_lower dem_vote_per if treat, lc(gs10) lp(dash)) ///
		(line y_hat dem_vote_per if !treat, lc(dknavy)) ///
		(line y_hat dem_vote_per if treat, lc(dknavy)), ///
		xline(50, lc(gs5)) ///
		legend(off) ///
		ti("`year'") ///
		xsc(r(25 75)) xlab(30(10)70) ///
		ysc(r(0 1)) ylab(0(.20)1) ///
		yti("Anti-Sanctuary Vote") xti("Dem Vote Share (%)") ///
		bgcolor(white) graphregion(color(white)) ///
		name("imm_voting_rd_plot_`yv'", replace)
	
	* Bring back the full dataset
	restore
}

* Put the plots together
graph combine imm_voting_rd_plot_2007_485 imm_voting_rd_plot_2012_366 ///
	imm_voting_rd_plot_2013_195 imm_voting_rd_plot_2017_342, ///
	graphregion(color(white))

